#include <stdio.h>
#include <math.h>
#define N 19
//找到指定点所在的区间
int find(double x[],double z)
{
	int i;
	int t = -1;
	for(i = 0;i < N-1;i++)
		if(z >= x[i] && z <= x[i+1])
			t = i;
	return t;
}
//函数值计算
double func(double x[],double y[],double h[],double m[],double z)
{
	int t;
	double s = 0;
	t = find(x,z);
	s = m[t]*powf((x[t+1]-z),3)/6/h[t] + m[t+1]*powf((z-x[t]),3)/6/h[t]+(y[t]-m[t]*h[t]*h[t]/6)*(x[t+1]-z)/h[t]+(y[t+1]-m[t+1]*h[t]*h[t]/6)*(z-x[t])/h[t];
	return s;
}
//一阶导数值计算
double gunc(double x[],double y[],double h[],double m[],double z)
{
	int t;
	double s = 0;
	t = find(x,z);
	s = -m[t]*powf((x[t+1]-z),2)/2/h[t] + m[t+1]*powf((z-x[t]),2)/2/h[t]+(y[t+1]-y[t])/h[t]-(m[t+1]-m[t])*h[t]/6;
	return s;	
}
//二阶导数值计算
double hunc(double x[],double y[],double h[],double m[],double z)
{
	int t;
	double s = 0;
	t = find(x,z);
	s = m[t]*(x[t+1]-z)/h[t] + m[t+1]*(z-x[t])/h[t];
	return s;	
}

int main()
{
	int n = N;
	double x[N] = {0.520,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.50,364.4,416.3,468,494,507,520};
	double y[N] = {5.288,9.4,13.84,20.20,24.90,28.44,31.10,35,36.9,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2};
	double z[5] = {2,30,130,350,515};	
	double y0 = 1.86548;
	double yn = -0.046115;
	double u[N];
	double v[N];
	double d[N];
	double h[N-1];
	double m[N];
	double t;
	int i;
	for(i = 0;i < n-1;i++)
		h[i] = x[i+1] - x[i];
	for(i = 0;i < n;i++)
	{
		if(i == 0)
		{
			v[0] = 1;
			d[0] = 6.0/h[0]*(y[1]-y[0]-y0);
			u[0] = 0;
		}
		else if(i == n-1)
		{
			v[n-1] = 0;
			d[n-1] = 6.0/h[n-2]*(yn-(y[n-1]-y[n-2]));
			u[n-1] = 1;			
		}
		else
		{
			v[i] = h[i]/(h[i-1]+h[i]);
			d[i] = 6*(y[i-1]/h[i-1]/(h[i-1]+h[i])+y[i+1]/h[i]/(h[i-1]+h[i])-y[i]/h[i-1]/h[i]);
			u[i] = h[i-1]/(h[i-1]+h[i]);
		}
	}
	for(i = 1;i < n;i++)
	{
		t = u[i]/2;
		u[i] -= t*v[i-1];
		d[i] -= t*d[i-1];		
	}
	m[n-1] = d[n-1]/2;
	for(i = n-2;i >= 0;i--)
		m[i] = (d[i]-v[i]*m[i+1])/2;
	//5个点的函数值：
	printf("5个点的函数值\n");
	for(i = 0;i < 5;i++)
		printf("%f %f\n",z[i],func(x,y,h,m,z[i]));
	//5个点的一阶导数值：
	printf("\n5个点的一阶导数值\n");
	for(i = 0;i < 5;i++)
		printf("%f %f\n",z[i],gunc(x,y,h,m,z[i]));
	//5个点的二阶导数值：
	printf("\n5个点的二阶导数值\n");	
	for(i = 0;i < 5;i++)
		printf("%f %f\n",z[i],hunc(x,y,h,m,z[i]));
	return 0;
}
